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Abstract 

When the Hamiltonian if of a system is represented by a finite matrix if, constructed 
from a discrete basis with overlap matrix B_, the matrix representation G(E) = 1/ {EB_ — if) 
of the resolvent G(E) = 1/(E — if) covers only one branch of G(E). We show how all 
branches can be specified by the phase of a complex unit of time = toe 1 ^. This permits 
ff to be constructed on a real basis; the only duty of the basis is to span the dynamical region 
of space, without regard for the particular asymptotic boundary conditions that pertain to 
the problem of interest. Specifically, we show that P cont G(E), where P con t projects onto the 
continuous spectrum of if, has the series representation 

(1 00 1 \ 

- + 2it/ff £ -X n (E^)4 1 2 1 (2^if) , 
L n=l 71 J 

where both the associated Laguerre polynomials, L^l^t^H), and the coefficients T n {Et$), 
satisfy 3-term recurrence relations. 



1 Introduction 



The theoretical treatment of a continuous stationary or quasistationary process requires 
some knowledge of the resolvent G(E) = 1/(E — H), where H is the Hamiltonian of the 
system at hand. There are numerous methods for evaluating G(E), at least approximately. 
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Frequently H is approximated by a finite-dimensional matrix I£. The straightforward sub- 
stitution of H_ for H leads to the approximation of G(E) by the matrix G(E) = 1/ (EB — H) 
where B_ is the overlap matrix for the discrete basis functions used in constructing H_ . In 
principle, the basis functions are only required to accurately span the (generally finite) spa- 
tial region where the dynamics take place. However, in practice, unless the wavefunction is 
matched to its asymptotic form outside the dynamical region, as is done for example in the 
R-matrix method, |l| the basis may also be burdened by the requirement that it incorporate 
the asymptotic boundary condition. This is evident from the absence of branch points in 
G_(E); the expression 1/(EB_— H_) has only poles (since the eigenvalue spectrum of H_ is 
discrete). Some of these poles (those in the bound-state region) are legitimate, while the 
others (those in the scattering-state region) simulate the branch cut of G(E) by a sequence 
of discrete points. Hence G(E) represents, at best, only one of the branches of the exact 
resolvent, and this particular branch is determined by the choice of basis. 

In this paper we derive a series representation of the resolvent which is such that the 
asymptotic boundary condition can be imposed independently of the basis, without matching 
the wavefunction; a particular branch of G(E) is specified by the phase of a complex unit 
of time. In a companion paper we resum this series to give an integral representation 
of G(E), which is perhaps simpler to use in numerical applications. Both the series and 
integral representations can be used even when the Hamiltonian is approximated by a real 
Hermitian matrix H_, constructed from a real basis. This flexibility may permit a real 
basis to be used more generally and efficiently than otherwise possible, and, in addition, 
it may prove to be useful in the development and practical realization of a general theory of 
continuum processes in which detailed knowledge of the asymptotic form of the wavefunction 
is unnecessary. Recently a reformulation of perturbation theory for multiphoton processes 
was given || which allows, at least in principle, rates and branching ratios to be calculated 
in terms of the flux through a large hypersphere, without input on the asymptotic form of 
the wavefunction. However, the practical efficiency of this reformulation depends on the 
availability of a tractable method for representing the resolvent without explicit reference to 
the asymptotic boundary condition. 

We pause for a moment to mention some of the difficulties that arise when a discrete 
basis is required to simulate the asymptotic behavior of the open channels of a system. For 
the sake of discussion we consider a specific but typical basis that is used to construct the 
Hamiltonian of a one-electron system such as a hydrogen atom; the basis functions are 

v /z 2m(2mr)' +1 e iKr P n (2zfi:r)Y im (x), 

where P n (x) is a polynomial of degree n, where Y; m (x) is a spherical harmonic, and where 
k is a parameter — the wavenumber of the basis. This basis can be readily generalized to 
a multiparticle system by specifying the orientation of the system in terms of Euler angles 
and by introducing products of radial functions, one in each of the interparticle distances 
(or combinations of them), perhaps associating a different wavenumber with each distance. 
A bound-state wavefunction satisfies a real damped-wave asymptotic boundary condition, 
which can be described by choosing a real basis with positive pure imaginary wavenumbers. 
Hence a real basis is well-suited to the description of closed channels and to the study 
of bound-state properties. 0] However, when G(E) is constructed from a real basis it has 
spurious poles on the "unitarity" branch cut of G(E) along the positive real energy axis, and 
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these poles give rise to spurious resonances with zero width in the scattering-state energy 
region. The wavefunction for a compound (metastable) state satisfies a complex outgoing- 
wave asymptotic boundary condition, which can be described by choosing a complex basis 
with wavenumbers lying in the upper-right quadrant of the wavenumber plane; [] since the 
factor e lKr behaves as a damped outgoing-wave, this basis is well suited to the description of 
both closed channels and cmtgoing-wave open channels, and can be used to treat half-collision 
process, e.g. autoionization || or photoionization |J. Furthermore, G(E) now has poles in 
the scattering-state region that lie in the upper half of the complex energy plane, and they 
simulate a branch cut along a line which is distinct from the real axis, so spurious resonances 
are only of minor significance. However, a price is paid for choosing the basis to be complex, 
i.e. H_ is nonHermitian, and therefore unitarity is only approximately satisfied. In fact, the 
norm of a state vector is a nonanalytic function of K since the bra (\&| depends on 

k* ; consequently, the norm does not converge as the basis size is increased. In addition, it is 
difficult to extract partial rates, i.e. branching ratios, since a partial rate is also nonanalytic 
in k. Finally, a complex basis cannot, in general, describe the real standing-wave boundary 
condition satisfied by the wavefunction for a full-collision process. Q 

While a real basis, with pure imaginary wavenumbers, can describe a standing-wave 
boundary condition, the branch of G(E) cannot be specified uniquely by real basis functions 
since they do not distinguish between ingoing- and outgoing-wave behavior. On the other 
hand, when a real basis is used in conjunction with the series and integral representations 
developed here, the only duty of the basis is to span the dynamical region of space. The 
only purpose of the purely decaying exponential factor e lKr (k is positive pure imaginary) is 
to restrict the range of the basis to the dynamical region. 

In order to impose asymptotic boundary conditions that are independent of the basis we 
start from the observation that all the particles in a system evolve according to a common 
time, and the asymptotic behaviour of the system's wavefunction emerges from the initial 
boundary conditions, at the start of the system's evolution, say time t — 0. To simplify our 
analysis we consider a system comprised of only one particle; the generalization to more than 
one particle is fairly straightforward and will be dealt with elsewhere. The temporal behavior 
is governed by the time-evolution operator U(t) = e~ tHt (we set ft = 1 throughout). Let us 
introduce the dimensionless variable r = t/t , where to is the unit of time that characterizes 
the time scale on which the motion of the system occurs. If Im E > we can represent 
G{E) by an integral along the positive r-axis: 

poo 

G(E) = -it / dT e i(toE)T e -i(toH)r_ ,y 

Jo 

We can analytically continue this representation of G(E) to a sector of the .E-plane that 
includes the negative real .E-axis if we first project out the bound states, i.e. if we remove the 

1 Choosing the basis wavenumbers to lie in the upper right quadrant of the wavenumber-plane is equivalent 
to rotating the particle coordinates into the lower right quadrant of the complex position-coordinate-plane. 
See e.g. §. 

2 An exception occurs when the scattering potential falls off at least as fast as an exponential potential. 
In this case the deviation of the wavefunction from a standing wave is inconsequential at large distances. A 
wavefunction that is a real standing wave over short distances has a finite number of nodes, and it can be 
constructed from a complex basis provided that polynomials P n (2iKr) of sufficiently high degree are included. 
At the same time the correct branch of G_(E) is guaranteed by the factor e lKr . See e.g. ||. 
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bound-state poles of G(E). Thus we consider P cont G(E)P cont , where P con t is the projection 
operator P con t = 1 — Pbd, with P bd the bound-state projection operator. We rotate the contour 
of integration from the positive real r-axis through a small angle into the lower-right quadrant 
of the complex r-plane, and thereby obtain a representation of PcontG(P)Pcont which is valid 
in a region including the negative real P-axis. We can represent G(E) throughout the lower 
half of the P-plane by an integral along the negative r-axis: 

G(E) = -it / dT e i(toE)r e -i(t H)r^ ^ 

Jo 

and again by considering P cont G(P)P cont , and rotating the integration contour through a 
small angle (now into the lower-left quadrant of the r-plane) we obtain a representation of 
Pcont G(E) P C ont which is valid in a sector of the P-plane that also includes the negative real E- 
axis. As we let E approach the positive real axis these two different representations yield two 
different branches of P cont G(P)P cont , and since both representations give the same "physical" 
branch of P con tG(P)P C ont on the negative real P-axis we can define a global representation 
of P cont G(E) P C ont on the first — the "physical" — sheet of the Riemann energy surface by 
using Eq. (JJ) for Im E > and Eq. (0) for Im E < 0. The point is: A particular branch 
of Pcont G(E) P C ont on the positive energy axis can be specified by the contour of the global 
integral representation. Our goal is to perform the integration over r in such a way that a 
signature of the contour is preserved. 

That signature is the angle of rotation of the contour, whose vestige is the phase of a 
complex unit of time = t$e 1 ^ '. The main result of our paper is that the resolvent, with the 
bound-states removed, can be expressed generally as the series 

(1 CO 1 \ 

- + 2tt (j> 2 H Y, -2n(P^)Pi 1 2i(2^P")J , (3) 

where L^l^t^H) is an associated Laguerre polynomial of degree n — 1 in the operator 2t ( j ) H, 
and where the coefficient Z n (Pt<^) is a number defined by the integral 

l n (a) = fVeW— V. (4) 

JO \T — lJ 

The energy, P, can take on any value on the first sheet of the Riemann surface cut along 
the line arg (P) = — 0. The whole Riemann surface can be covered by rotating the cut, i.e. 
by varying the phase 0. The magnitude, to, of the unit of time, t,p, is a parameter whose 
value is to some extent, but not entirely, arbitrary. It should not differ greatly from the 
characteristic unit of time for the motion of the system, i.e. from the time it takes for the 
wavefunction of the system to change appreciably, for otherwise convergence of the series 
would be difficult to attain. In particular, to should not exceed the characteristic duration 
of the process under study, since a measurable process cannot take place in a time less than 
it takes for the wavefunction to evolve appreciably. For example, the characteristic duration 
of an elastic collision process at energy E is of order 1/P so that in this case we require 
Eto < 1. For asymptotically large values of E the time-scale is infinitesimally short, i.e. 
t$ ~ 0, and only the first term on the right side of Eq. (|j) contributes, giving G(E) ~ 1/P 
for all branches. Apart from this 1/P term, the dependence of G(E) on P is contained in 
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the coefficients T n {Et^), and a particular branch of G(E) is fixed by these coefficients after 
specifying the phase 0. The physical branch of G(E), which has outgoing- wave behavior, is 
specified by choosing to be in the range < < tt, and the other, "unphysical" , branch, 
which has ingoing-wave behavior, is specified by choosing to be in the range > > — ir. 
In principle, the right side of Eq. (|3]) should not vary as varies over its allowed range, 
but in practice, when the series is truncated and the summation not fully converged, it does 
vary — the less so, the better the convergence. The presence of P CO nt in Eq. (0) is perhaps 
not surprising since we were obliged to project out the bound states before writing down 
a global integral representation of the resolvent. [Actually, since G(E) acts, in general, on 
a wavepacket, we do not need to multiply the Hamiltonian by P CO nt; rather, we need only 
let P CO nt act on the wavepacket.] When the Hamiltonian contains a potential that has an 
attractive Coulomb tail, an infinite number of bound states accumulate at threshold, and 
it is neither necessary nor desirable to remove those bound states just below threshold; we 
elaborate further on this in sections 2 and 3. 

Note that the Laguerre polynomials satisfy a 3-term recurrence relation, and therefore 
they can be calculated recursively with only a minimal number of multiplications of H_. In 
addition, as shown in Appendix C, the coefficients T n {Et$) can be calculated using a 3-term 
recurrence relation. Furthermore, since the evolution operator e~ l <t> H acts for only one unit 
of time, it can be calculated using a Pade approximant — see Appendix C. 

In the derivation of Eq. (|3|) which we give in section 4, we initially restrict the phase of 
tfy to the range — n/2 < < ir/2. However, the series can be analytically continued outside 
this range. The preliminary restriction on the phase is related to a restriction on the contour 
of the global integral representation, as we now explain. For the sake of discussion, assume 
E to have a physically realizeable value, one appropriate to a collision process, i.e. E real 
and positive. The integrals of Eqs. (|l|) and (0) differ only by their contours. These contours 
can be rotated from the positive and negative real t-axes into the upper half of the t-plane, 
but they must remain distinct since the integrals represent two different branches of G(E). 
Hence if we choose one contour to lie along the ray arg (t) = in the upper-right quadrant, 
with therefore restricted to < < 7r/2, we should choose the other contour to lie along 
the ray arg (t) = n + in the upper-left quadrant, with restricted to — 7r/2 < < 0; 
therein lies our preliminary restriction on the phase. 

To gain an intuitive understanding of the series, we give now a brief heuristic derivation 
of Eq. (|3p. We start by formally expressing G(E) as a Taylor series in H/E: 



This series is an asymptotic series in powers of H/E. It is applicable for large E, i.e. it is 
applicable on a short time-scale, but not in all sectors of the complex P-plane. In fact, it is 
not applicable when E is real and positive; for if H is also positive, i.e. if H acts on those of its 
eigenvectors belonging to the continuous eigenvalue spectrum — whose threshold we choose 
to be zero — all of the terms (H/E) n are positive and the sum steadily diverges as more 
terms are included. Hence, for E real and positive we are motivated to delete continuous 
eigenvalues much larger than E by multiplying the right side of Eq. (0) by a cutoff factor. A 
natural choice for the cutoff factor is e~ l,t > R ', provided the following restrictions are imposed 




(5) 
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on tf. In order for t~ l 4> H to cut off the sum we require that Re > 0, i.e. \<f>\ < tt/2, 
but since we do not need (or want) to delete continuous eigenvalues much smaller than E 
we require that Et < 1; these are the same restrictions on mentioned above. (We note 
again that the restriction on can be lifted by analytic continuation.) We can regather the 
terms in the sum on the right side of Eq. (|5|), now modified by the cutoff factor, to yield 
an expansion in the basis functions xe~ x ^ 2 L^\x), where x = 2t^H and where the L^\x) 
are associated Laguerre polynomials. These basis functions satisfy the end-point boundary 
conditions for the modified sum, i.e. they vanish linearly at x = and exponentially at 
x = oo. Furthermore, these basis functions form a complete set over the interval < x < oo, 
and they satisfy the orthogonality relation 

f° dx x v e- x Lt\x)L^\x) = ^±^±il j mn , ( 6 ) 
Jo 1 [n + 1) 

Here x is understood to be a continuous variable, and consequently we must restrict the 
operator H to act only on those of its eigenvectors belonging to the continuous eigenvalue 
spectrum. In other words, we must exclude eigenvectors belonging to the discrete eigenvalue 
spectrum. Choosing the upper index v of the Laguerre polynomials to be 1, it follows that 

P '""" " Pconte"*^ ( i + it^H) f) -l^EQL^l^H)) , (7) 



E-H ^ \E " r r " 
where the coefficients, determined using Eq. (Bl), are 



oo 



IjEt..): -2,/ dHe-^l^^-^jL^^H), (8) 



o 



where here H is understood to be a nonnegative continuous variable of integration. Thereby 
we arrive at Eq. ([|). While not obvious, the coefficients defined by Eq. (§) are the same as 
those defined by Eq. (|j); we prove this in Appendix A. Although Eq. (|7|) was just derived 
for E real and positive, the results can be analytically continued to complex E. 

It is time to remark on other polynomial expansions of the resolvent. If the Hamiltonian 
is represented by a finite matrix, whose minimum and maximum eigenvalues are E min and 
-Etaax; respectively, one can, by fiat, expand the matrix representation of the resolvent in 
any set of orthogonal polynomials that form a complete set over a finite interval rescaled to 
the interval [£ min ,£ max ]. The most popular of these expansions uses Chebyshev or Faber 
polynomials [^|, |TIJ since the evolution operator U (t) = e~ lHt can be expanded very effectively 
over a finite time interval in terms of these polynomials. |T0|, jlTf The Chebyshev expansion 



has been remarkably successful in obtaining properties of bound and compound states. JT^ ] 
However, the Chebyshev and Faber expansions do not discriminate between the discrete 
and continuum parts of the spectrum of H, and notwithstanding some numerical evidence 
that these expansions can yield reasonably accurate estimates of scattering amplitudes in 
some simple cases, [|, 0] it is unclear at this stage how well-suited are these expansions 
for the treatment of general collision processes, particularly when the Coulomb tail plays an 
important role. The long-time behavior of the temporal correlation amplitude (see below) is 
determined by the spectral density at the threshold of the continuous spectrum. The success 
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of the Chebyshev expansion in dealing with discrete (bound and compound) states can be 
understood by observing that each of these states decays exponentially with increasing time 
in a wide sector of the complex time-plane; hence the effective time-interval over which 
U(t) must be expanded is relatively short. However, in a continuum process, wavepacket 
spreading occurs, and its asymptotic temporal behavior is 1/t 3 / 2 (for short-range potentials), 
a rather slow fall-off with t. 

Our primary interest is in stationary processes in the continuum, for which E is real and 
positive. We focus in the present paper on the general theory underlying the series represen- 
tation of the resolvent. As noted already, for simplicity we consider a system comprised of 
only one particle. Our analysis is largely informal (a mathematically rigorous treatment is, 
anyway, beyond our abilities). In a companion paper we derive the integral representa- 
tion, and we give numerical illustrations for the examples of photoionization of a hydrogen 
atom and s-wave scattering from a 1/(1 + r) 4 potential. We plan to report on an application 
to a multiparticle system in the future. 

We frame our paper within the context of the inclusive rate at which a continuous sta- 
tionary or quasistationary process occurs. This rate can be expressed in terms of a matrix 
element of G(E) of the form {ip\G(E)\ip) where is a localized wavepacket. This 
matrix element can in turn be expressed, according to Eq. (|l]), as — i / °° dt e lEt C(t) where 
C(t) is the temporal correlation amplitude: 

C(t) = (^(t)), (9) 

with \ip(t)) = U(t)\ip). In the next section we explore the analytic properties of the correla- 
tion amplitude in the complex-time-plane. We find that the singularities of C(t) lie in the 
upper half of the t-plane. While we can expand U(t), and thereby C(t), in powers of t, the 
power series for C (t) has only a small (maybe infinitesimal) radius of convergence. In section 
3 we transform variables, from t to a variable u, to obtain a power series which is more useful. 
Thus we make a conformal transformation which maps a singularity-free region of the t-plane 

- namely, the half of the complex t-plane that lies below the line Im t/Re t = tan0 - 
into the unit circle in the w-plane. We can express U(t) and C(t) as power series in u, with 
coefficients c n that are rather simply related to the coefficients of the power series in t. We 
analyse the properties, in particular, the large-n behavior, of the coefficients c n , though we 
defer some of the analysis to Appendix D. In section 4 we perform the integration over t and 
obtain the series representation, Eq. (|j), of G(E). We analyse the convergence properties of 
this series, and also show that the higher terms in the series can be resummed as another 
series that converges rapidly. In Appendix A we prove the equivalence of the two forms of 
X n (a) introduced above. In Appendix B we use the series representation of the resolvent to 
reproduce some known formal results, for example, threshold laws. In Appendix C we de- 
scribe some algorithms that are useful for the implementation of the series, and in Appendix 
D we analyse the large-n behavior of the coefficients c n when the potential has a Coulomb tail. 

2 Analytic Properties of the Correlation Amplitude 

Recall that is a normalizeable, localized wavepacket which evolves in time t as 
= U(t)\ip), where U(t) = e~ iHt is the t ime-evolution operator and where H is the 
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Hamiltonian of a system that is comprised of only one particle, of mass /j, say. We assume 
that the wavepacket \ip) includes a continuous superposition of scattering eigenstates of H, 
and is not merely a superposition of discrete (bound- and compound-) eigenstates. If the 
potential has an attractive long-range Coulomb tail, a localized wavepacket has, in general, 
a nonzero overlap with an infinite number of bound-eigenstates of H. Indeed, if L is the 
characteristic linear dimension of \ip) in position space, Rydberg bound states with energy 
eigenvalues greater than or of the order of —n/L are indistinguishable in the composition of 
the wavepacket from scattering eigenstates with energy eigenvalues less than or of the order 
of n/L. As a consequence, high Rydberg bound states, through their contribution to the 
correlation function, play an important role in continuum processes, a feature we observe 
below. The case where the potential has a Coulomb tail is usually an exception requiring 
special treatment. Leaving aside this case for the moment, the continuum portion of \^>{t)) 
spreads linearly in time, and therefore occupies a volume (in 3-dimensional position space) 
that is proportional to t 3 . Since U(t) is unitary, (if)(t)\ip(t)) is conserved in time, and so 
the continuum portion of the wavepacket attenuates at each point in space as t~ 3 / 2 . Conse- 
quently, the correlation amplitude C(t) = (ip\ip{t)) contains a component which vanishes as 
£-3/2 f or f oo. Therefore C(t) has a branch point at infinity, and since this branch point 
is of order two there must be another branch point, joined to the branch point at infinity by 
a cut in the complex t-plane. 

As an example, consider a free particle (of mass /x) whose position is initially described 
by the Gaussian wavepacket 



The Hamiltonian governing the evolution of the wavepacket is the kinetic energy operator 
-(l/2/i)V 2 ; we find that 



where we recall that to characterizes the time scale for the evolution of the wavepacket 
and is defined in the present example as to — /V K o- Note that C(t) has branch points at 
t = 2Uq and t — oo. At time t — the wavepacket has a characteristic width in position 
space of 1/k,q, but since the wavepacket has a momentum distribution of width kq its spatial 
spread after time t is K t//i, and this spread exceeds the original width of the wavepacket 
when t is comaparable to t . In other words, the singularity at 2it is a signature of the 
time at which the wavepacket becomes significantly deformed by spreading. The wavepacket 
can evolve either forwards or backwards in time to the single point t — oo, and the result 
depends on the arrow of time. Hence the wavepacket is double-valued at t — oo, a property 
that is encompassed by the branch point at t — oo. 

A general (nonGaussian) wavepacket has a momentum distribution whose very-high- 
momentum components are appreciable, and the spatial tail of the wavepacket deforms 
the moment it begins to evolve freely. Even the tail of a Gaussian wavepacket deforms 
instantly if it evolves under the influence of a potential that can transfer a large momentum 
to the particle. Hence, in general to — 0, i.e. the correlation amplitude has a branch 
point singularity at the origin. To be more concrete, suppose that H has bound-state and 
scattering-state eigenvectors |xbd,n) and |xk), respectively, with real energy eigenvalues E hdtn 




(10) 




(11) 
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and Ek = k 2 /2fi, respectively. We can express \ip) as the superposition 

|V>> = E^bd,n|Xbd,n) + / d 3 k ^(k)|xk), (12) 

n J 

where the eigenvectors are normalized so that V'bd.n — (Xbd^lV') and ip(k) = (%k|V')- This 
superposition evolves in time as 

m)) = E^bd,n|Xbd,n)e- 1 ^ + / rf 3 A; V(k)e-^*| Xk ), (13) 

n J 

and the correlation amplitude is 

C(t) = C hA (t) + C amt (t), (14) 

where 

C bd (0 = El^d,n| 2 e-^ d -', (15) 

n 

C cont (t) = /rf 3 A;|^(k)| 2 e-^. (16) 

We can let t move into the lower-half complex t-plane; the exponential e~ lEkt decays with 
increasing k and both Ct>d(£) and C cont (t) are well-defined. However, e~ tEhd - nt explodes as 
t — > oo in the lower-half of the t-plane; hence Cbd(t) has an essential singularity at t = oo 
and is unbounded. Furthermore, if we allow t to move into the upper-half complex t-plane, 
e ~*E k t eX pi oc i es with increasing k, as the Gaussian exp(i?fc Im t). Therefore, unless |?/>(k)| 2 
decreases more rapidly than this, C cont (t) is formally undefined in the upper-half of the t- 
plane. Suppose for the moment that |^(k)| were to decrease with increasing k as e~ Ekt °. In 
this case C con t(£) would be formally defined in the region Im t < 2t . However, by rotating the 
contour of /c-integration through an angle 6 into the first octant of the lower-right quadrant 
of the A;-plane, assuming that 1^(^)1 is free of singularities in this octant, we can analytically 
continue C cont (t) throughout any finite region of the sector < arg(t) < 20 in the upper- 
right quadrant of the t-plane, excluding the section of the positive imaginary t-axis above 
2it since both e~ EkU) and e~ tEkt are undamped oscillatory functions of k when 6 = 7r/4 
and t is pure imaginary. Therefore C cont (t), and hence C(t), are analytic in both the lower- 
and right-half t-planes, but they have branch points at 2Uq, and branch cuts extending from 
2ito to infinity. In general, IVK^OI decreases as a power of \ jk with increasing k, less rapidly 
than a Gaussian, and we may expect C(t) to have a branch point at the origin, i.e. a branch 
point at 2it with t ~" *■ 0. Furthermore, (^(k)] generally has singularities at finite points 
in the complex /c-plane (in contrast, e~ Ekt ° has an essential singularity at k — oo); but if 
there were no singularities in the lower-right quadrant of the /c-plane, we could rotate the 
contour of /c-integration so that it runs along the right edge of the negative imaginary axis, 
and we could subsequently move t from the positive real axis through the upper-half t-plane 
to the upper edge of the negative real axis. In this case C(t) would be defined everwhere in 
the finite complex t-plane, except possibly for a branch point at the origin, which we cannot 
exclude unless we can move t continuously along a closed loop around the origin (without 
discontinuously moving the contour of /c-integration). Hereafter we formally define t such 
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that the singularity of C(t) nearest to the origin is located at 2Uq, and we draw a cut along 
the positive imaginary axis from 2zt to oo. Of course, if, by this definition, i were to vanish, 
it could not serve as a unit of time; however, this is not of concern in practice, for in practice 
to may be small but it does not vanish, a matter we return to at the end of this section. 

To determine the general asymptotic behavior of C(t) for t ~ oo let us transform the 
integration variable k to k/ \/t on the right side of Eq. (|T6|). Assuming the potential is short- 
range, we can, as a first approximation, replace \ip(k/\/t)\ 2 by \ip(0) | 2 for t ~ oo. However, 
if the wavepacket carries orbital angular momentum |?/>(0)| 2 may vanish; small values of k 
are inhibited by the centrifugal barrier. If I is the smallest angular momentum quantum 
number present in the wavepacket, |?/>(k)| 2 vanishes as k 21 for k ~ 0. Therefore we write 

|^(k)| 2 ^hMk)| 2 . (17) 

It follows that 

C C ont(t) ~ e-^ +3 ^ 4 (2vr) 3 / 2 (2Z + 1)!!|^(0)| 2 Wtf+ 3/2 , t „ ^ (1 8 ) 
where we used 

Jo 2(2a) m V a 

Note that the large-t behavior of C cont (i) is determined by small k, i.e. by the continuum 
eigenvalues of H that are close to threshold. If the potential has a Coulomb tail |^(0)| 2 
is infinite, and we must factor out the divergence arising from the normalization constant. 
Thus, if the Coulomb tail is —Ze 2 /r at large radial distances r, we remove an offending 
factor from |?/.>(k)| 2 by writing (for all I) 




where 7 = Z/(aok), with ag = l//xe 2 . After making the transformation k — > k/\/i, we have 
7 ~ (Z/aok)\/i and the prefactor on the right side of Eq. (|19D either vanishes exponentially 
(if Z < 0) or becomes the divergent function {2nZ / aok)\fi (if Z > 0) when we let t increase 
to infinity in any sector excluding the negative real axis. Hence, when an attractive [Z > 0) 
Coulomb tail is present, C cont (t) falls off only as \jt: 

Ccont(i) ~ -8z7r 2 Z/i|^(0)| 2 /(a t), t~oo (Z>0). (20) 

Note, however, that a Coulomb tail gives rise to an essential singularity at t — 00. (If we 
were to let t increase to infinity in any sector excluding the negative real axis on the second 
sheet of the Riemann t-surface, C con t(t) would vanish exponentially if Z > 0.) A branch 
point at t = 00 remains, but it does not dominate the asymptotic behavior. Futhermore, if 
Z > 0, Rydberg states converging to threshold from below cannot be distinguished in the 
wavepacket from continuum eigenstates converging to threshold from above, and therefore 
we have Chd(t) ~ C cont (t) for t approaching 00 in the upper-half t-plane (the half-plane 
in which lower-lying bound states decay exponentially). It follows that, in the upper-half 
t-plane, 

C(t) 16m 2 Zti\i>(0)\ 2 /(a t), t ~ 00 (Z > 0). (21) 
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Figure 1: The contour C runs along the upper edge of the real axis of the physical sheet of 
the Riemann energy surface cut (the zig-zag line) along the positive real axis (the "unitarity" 
cut). Bound state poles, indicated by x, lie on the negative real axis of the physical sheet, 
while resonance poles, indicated by *, lie on the unphysical sheet reached by crossing the 
cut. The resonance poles are distributed symmetrically about the cut. The contour C is the 
result of bending C around the negative imaginary axis; the left half of C lies on the physical 
sheet, while the right half lies on the unphysical sheet. 



In other words, we gain a factor of 2 in the asymptotic form of C(t) through the contribution 



of high Rydberg states. As seen in Eq. ( 122 ) of Appendix B, this factor of 2 enters the rate 
for a continuum process near threshold. When the potential has a Coulomb tail, the sum 
(over bound states) on the right side of Eq. (|T5| ) is convergent, but not uniformly, for t in the 
upper-half t-plane, and we cannot interchange the limit t — > oo and the sum. In contrast, 
when the potential is short-range C^{t) ~ for t ~ oo in the upper-half t-plane. 

To gain further insight into the behavior of C(t) at large t it is useful to express the 
time-evolution operator U(t) in terms of the resolvent G(E). We have |]13 



U(t) = 7T-. I dEe- iEt G(E), (22) 

where, assuming that t is real and positive, the contour C runs along the upper edge of the 
real E-axis from oo to — oo; see Fig. 1. The resolvent has branch points at E — and E = oo, 
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and is defined on a two-sheeted Riemann energy surface; the "unitarity" branch cut is drawn 
along the positive real E-axis. In addition to branch points, G(E) has bound-state poles at 
points -Ebd,n on the negative real energy axis of one sheet — the physical sheet — and has 
resonance poles at points E res n and E* esn in the lower and upper half-planes, respectively, 
of the other sheet — the unphysical sheet. Let us bend the contour C around the branch 
point at E — 0, into the lower-half i?-plane, so that the new contour, C, wraps around the 



negative imaginary axis. |13fl As we distort C it sweeps over the bound-state poles (on the 
physical sheet) and also over those resonance poles dispersed in the lower-right quadrant on 
the unphysical sheet of the energy-plane. It follows that 

\m) = 1^— J c dE e~ iEt G(E) (23) 

= E^bd,n|Xbd,n)e- iBbd -' + £ ^ reSl „|Xre Sj n}e- iiW + |^g(t)>, (24) 
n n 

where |Xres,n) is an eigenvector of H satisfying out going- wave boundary conditions corre- 
sponding to a compound (resonance) stated with a complex energy E Tes ^ n whose real part is 
positive, where ip res ,n = (Xres.nlV'); an d where |^bg(^)) describes the continuum background: 

\4> hg (t)) = ^-J c ,dEe-* Et G(EM). (25) 

As expected from our earlier discussion, if the potential is short-range the vector |^b g (^)) 
attenuates as t -3 / 2 with increasing t; a general proof, due to Zumino, is presented in Ref. fT3|| . 
In the case where the potential has a Coulomb tail, the behavior of \ipbg{t)) at large t has 
been analyzed by Dollard. [fyj Note that once C has been deformed to C we can analytically 
continue the right side of Eq. (p5|), and hence C(t), from the real positive t-axis to the entire 
right half of the complex t-plane. 

When t is real and negative, U(t) has an integral representation similar to the right side 
of Eq. (^) but with a contour C running along the lower edge of the real i?-axis and in the 
direction opposite to C. Assuming that the Hamiltonian is invariant under time-reversal, the 
two representations are related through 

G(E*) = KG(E)K\ (26) 

where K is the antiunitary time-reversal operator. |1| Hence, if represents the time- 
reverse of the wavepacket i.e. if = K\ip), we have, with t real and negative, 

\m) = ^~ j/E e- iEt G(E)$) (27) 

= — / dE e- iEt KG(E*)K^\$) (28) 
2iri Jc 

= K-^— [ dE e iEt G{E*)\ip) (29) 
2m Jc 

= K\iP(-t)), (30) 



3 When Ek is complex the continuum eigenfunctions (x|xk) explode exponentially in position space, as 
exp(ik • x), for |x| ~ oo. Consequently, certain integrals over x — for example, ("0|Xrcs,n) — are formally 
undefined. However, such integrals may be defined through analytic continuation, e.g. consider / °° dr e ar , 
which is — 1/a for all a ^ 0. 
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where in the second step we used Eq. fl2"S| ) and in the third step we noted that K^K = 1 
and that K complex-conjugates c-numbers. If on the right side of Eq. fl27|) we deform the 
integration contour into one wrapped around the positive imaginary E'-axis, we obtain an 
expression for \ip{t)) that is similar to the right-side of Eq. ( p4|) but with all terms time- 
reversed. The time-reversed resonance terms correspond to the conjugate poles of G(E) in 
the upper-right quadrant of the .E-plane on the unphysical sheet. Let us introduce the new 
correlation amplitude 

C(t) = (V#(*)>. (31) 
As long as t is real and negative we can use Eq. ( |3"0D to write 

C(t) = @\ [K\i/,(-t))] (32) 

= (iP(-t)\[K^)] (33) 

= m-t)\\^) (34) 

= [C(-t)Y, (35) 



where in the second step we noted [15 that, since K is antilinear, (b\ (K\a)) = (a\ (K'\b)) 



for any two kets \a) and \b). After deforming the contour C, we can analytically continue 
C(t) into the entire left half of the complex t-plane, and since both [C^t*)]* and C(—t) are 
also analytic functions of t in this region we can generalize Eq. ( [3~5|) to 

[C(t*)}* = C{-t) (36) 

for t anywhere in the left half of the complex t-plane. When t lies on the negative imaginary 
axis, C{t) is real, and therefore [<?(-**)]* = C{t) = [C(-t*))*. It follows that C{t) is the 
analytic continuation of C(t) into the left half of the complex t-plane. If the point t is moved 
over the branch cut along the positive imaginary axis, from the left edge of the cut to the 
right edge, on the same sheet, the correlation amplitude jumps discontinuously from C(t) to 
C(t); in other words, at the cut C(t) and C(t) are different branches of the same multivalued 
correlation amplitude. 

To conclude: In general, C(t) is a nonsingular function of t in the right-half of the finite 
complex t-plane, but has a branch point on the positive imaginary axis at t = 2Uq. We can 
expand C(t) in the right-half of the t-plane as — c.f. Eq. (|HD — 

C(t) = CUt) + CrcsW + C bg (t) (37) 
= E l^bd,n| 2 e-^ d -* + E \^,n?z~ iEies ' nt + (V#bg(*)>, (38) 

n n 

where Cb g (t) falls off as t~ 3 ^ 2 or, if the potential has an attractive Coulomb tail, as 1/t for 
t ~ oo. If the potential has a Coulomb tail C^ g (t) has an essential singularity at t — oo. Due 
to the exponential factor(s), CUt) an( i C Tes (t) also have essential singularities at t = oo. In 
increases to infinity in the lower half of the the t-plane Cbd(^) explodes exponentially, 
and C ies (t) exhibits similar behavior for t within some other sector of the t-plane. However, 
Cbd(t) is bounded in the upper-half of the t-plane. Also, C res (t) is bounded in the lower- 
right quadrant of the t-plane since Re E res ^ n > and Im E res ^ n < 0. In fact, C Tes (t) is 
bounded in an even wider region since G(E) has no poles within a sector of the E-plane, say 
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—S < arg (E) < S , which contains the positive real energy axis; therefore Im E TCS ^ n t < 0, 
and C res (t) is bounded, throughout the sector < arg (t) < E of the t-plane. Hence C(t) 
is bounded throughout the sector < arg (t) < S . Similarly, C(t) is nonsingular in the 
left-half of the finite complex t-plane, C res (t) is bounded throughout the lower- left quadrant 
of the t-plane, and beyond, and C(t) is bounded throughout the sector tt — S < arg (t) < ir. 
We can generate a power series for C(t) by expanding U(t) = e~ tHt in powers of t; we have 

m) = g tlEpl^ (39 ) 

and therefore introducing the dimensionless variable r = t/t , assuming that t ^ 0, we 
obtain 

oo m 

c{t) = Y,m-iHt,r\^)— (40) 

m=0 m - 

The radius of convergence of this expansion of C(t) in powers of r is 2 since the singularity of 
C(t) that is nearest to the origin is, by definition, located at 2Uq. Note that we just assumed 
t 7^ 0, yet earlier we remarked that we can have t = 0. While in principle it is generally 
true that t = 0, in practice t ^ 0, as we now explain. The time that it takes the wavepacket 
to deform significantly from its initial form is governed by the highest speed in the velocity 
distribution of the wavepacket as it evolves in the presence of the potential. Now the highest 
speed in the velocity distribution is effectively determined by the largest eigenvalue, E max 
say, of the smallest matrix if that can accurately represent H in the expression e~ %m \^)). 
Hence, in practice the characteristic value of to is 1/E max , which may be very small but is 
nonzero. 

3 Conformal Transformation 

To calculate the rate for some continuous stationary or quasistationary process from a 
correlation amplitude we need to know the correlation amplitude for all t on either the 
positive or negative real axis. However, the power series in r is useful only for < |t| < 2t . 
Since C(t) is singular on the positive imaginary axis, it is expedient to divide the t-plane into 
two half-planes separated by a line through the origin, and to conformally map the lower of 
these half-planes (i.e. the one free of singularities) into the unit circle; see Fig. 2. Thus we 
change variables from t to 

u = —. (41) 

t - te l< P 

The mapping u depends on both t and the complex unit of time 

H = e l %, (42) 

and a particular branch of C(t) can be specified by 0, or, rather, a range of values of 0. As (j> 
varies over the range — n/2 < <fi < tc/2 the boundary line Im t/Re t = tan0 rotates through 
one revolution, and we remain on the same branch of C(t). To pass to another branch of 
C(t) we must allow <fi to move out of this range. We can do this by analytic continuation in 
the variable t^, to the left-half of the t^-plane, i.e. to values of in the range n/2 < |0| < ir. 
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a Imx 




Figure 2: The correlation amplitude C(t) has a branch-point singularity on the positive 
imaginary axis in the r = t/t plane, and another one at infinity. The unit of time to 
is defined so that the singularity nearest the origin is at r = 2i. We have drawn a cut 
extending up the positive imaginary axis from 2i. One-half (the hatched section) of the 
r-plane is conformally mapped onto a unit circle. 
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Figure 3: The unit circle in the w-plane, where u = (r + ie l ^)/(r — ie 1 ®). If the phase 
is positive (negative), the real positive (negative) axis in the r-plane is mapped onto the 
trajectory indicated by the broken line in the upper (lower) semicircle of the w-plane. This 
trajectory is furthest from the real w-axis at the point u or — u where u = i cos 0/(1 + 
sin |0|). The correlation amplitude has a singularity (indicated by •) at u = 1 and another 
singularity outside the unit circle. 



Note that the boundary line Im t/Ke t = tan0 is mapped onto the circumference of the 
unit circle in the u-plane, and the real positive (negative) t-axis is mapped onto the broken 
line shown in the upper (lower) semicircle of Fig. 3 if is positive (negative). The points 
t = 0, t = t$, t = oo, t = —ittj,, and t = it^ are mapped onto the points u = —1, u = i, 
u = 1, u = 0, and u = oo, respectively. There are no singularities inside the unit circle, 
and only one on the circumference, at u — 1 (corresponding to the singularity at t = oo). 
A singularity on the positive imaginary axis of the t-plane is mapped onto a point outside 
the unit circle in the lower (upper) half of the u-plane if is positive (negative), and moves 
onto the real w-axis as vanishes. 

We now express C(t) as a power series in u. Substituting for r in the power series on the 
right side of Eq. (|^) using 

T = - ie .*(I±^, (43) 



u 
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and noting that 



1 + tt 
1 — u 



m oo 



{mm 



n=0 



nl 



where P n {Tn) is the polynomial^ 

P n (m) 

we obtain the new power series: 



d" 



1+u 



du n V 1 — u 



u=0 



c(t) = £ c^k, 



n=0 



where 



and where, with z the operator —t^H, 



CJz) 



E 

m=0 



-2 ' 



n\m\ 



(44) 

(45) 

(46) 
(47) 

(48) 



Since u = — 1 when t — 0, and since C(0) = ('01V'}; we have 

oo 



(49) 



n=0 



which provides a useful check on the accuracy of the coefficients c n (t^) in a practical appli- 
cation. Since C(t) is nonsingular everywhere inside the unit circle in the u-plane, the power 
series in u, on the right side of Eq. (|46|), converges for all \u\ < 1. The point u = 1 requires 
special consideration, as discussed below. 

The series for CJz) can be expressd in closed form. This is obvious when n = since 
Pjm) = 1 and we have 

C (z) = e z . 



(50) 



For n > 1, we use Eqs. 



and 



to write: 



CJz) 



1 j d n 
n\ ] du n 



00 z m n + u ^ m 



m! V 1 — u 



.m=0 



w=0 



n! 



2z -S- 

-e 



(51) 



M=0 



Incidentally, it can be shown that 



(jTl ~\~ Tl — 1)! 

in (to) = — T -i— 2Fi(-n, -m, -m - n + 1;-1), n. > 1, 

(to — 1)! 

with Po( m ) = 1- The P n (m) are Krawtchouk polynomials, which for n > 1 satisfy the recursion relation 

P„(x) = 2a;P„_ 1 ( a + (n - l)(n - 2)P n _ 2 (z), 

with P (x) = 1- 
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where in the second step we noted 



l + u ' 



2z{ 



i-uj = e e u 



Now the generating function of the ordinary Laguerre polynomials, L m (x), is [16 



(1 - u)- l e~<^) = V L m (x)« m , (52) 



m=0 



and it follows that 



C„(z) = e z [L n (-2^) - L„_x(-2z)]. (53) 

Since L n (x)—L n _i(x) = — (x/njL^^x), where L^_ 1 (a;) is an associated Laguerre polynomial 
of degree n — 1, which may be recast as a confluent hypergeometric function n 1 F 1 (l — n, 2, x), 
we arrive at 

C n (z) = (2z) 1 F 1 (1 - n, 2, -2z)e z , n > 1, (54) 

which is the desired result. 

It is instructive to write C n (z) = <5n(<2)e z , where Q n (z) is a polynomial of degree n in z, 
so that 

Cn{U) = (i^lQni-UHMi-iQ), (55) 

where \ijj{—it^)) is the wavepacket that has evolved for the complex time —it^ from the 
wavepacket at t — 0; since the point i = — it^ is further from the singularity at 2zt than 
the point t = 0, the influence of the singularity on \if)(—it<f,)) is weaker than on The 
point £ = — it<f,, or, equivalently, r = —ie 1 ^, corresponds to the origin of the w-plane, i.e. the 
point about which a power series in u is developed. Hence the conformal transformation 
permits the analytic continuation of C (t) from a power series in r to a power series in u via a 
"connection" point r = —ie % ^ which lies within the circle of convergence of the power series 
in r, but on the side of this circle furthest from the singularity at 2Uq. 

We now explore some properties of the expansion coefficients (^(t^). It is evident from 
Eqs. (f47|), (|50D, and (|54]) that since H is Hermitian, and since C n (z) is a real function when 
z is real, 

[c n (t;)]* = c n {U). (56) 

Furthermore, since t*± = i_<^ we have 

C n (t-4,) = [Cn(W. (57) 

We know from the previous section that (if)\e~ t,l ' H \ip) is analytic in the t^-plane except on 
the negative imaginary axis, where = ±7r. Hence Cq(^) is analytic everywhere in the 
finite t^-plane cut along the negative imaginary axis. Since Q n (—t^H) is a polynomial, and 
therefore an analytic function of t^H, we infer that Cn(t^) is also analytic everywhere in the 
finite cut t^-plane, where n is any nonnegative integer. 

Although the expansion coefficients c n (t$) are analytic throughout the range \<p\ < tt, if 
we allow to move out of the range |0| <n/2 the boundary line Im i/Re t = tan0 crosses 
the cut that we have drawn along the positive imaginary t-axis. Hence ii tt/2 < \<p\ < ir this 
cut is mapped into the unit circle and the expansion of C(t) in powers of u, i.e. Eq. 
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longer converges for all u within the unit circle. Nevertheless, provided that u is sufficiently 
small, i.e. provided that t is sufficiently close to — it^, the expansion of C(t) in powers of 
u converges even for tt/2 < \(j)\ < tt. As increases from to 7r_ — where in general a_ 
and a + , respectively, are numbers just below and just above the number a — the point —it^ 
moves from the negative imaginary axis, into the right-half of the t^-plane, and onto right 
edge of the positive imaginary axis. Thus, there is always a region in the right-half of the 
t-plane within which the expansion of C(t) in powers of u converges when < < tt; but 
as approaches 7r_ this region shrinks to the point i on the cut in the t-plane. Similarly, 
as decreases from to — 7r_, the point — it^ moves from the negative imaginary axis, into 
the left-half of the t^-plane, and onto left edge of the positive imaginary axis. Thus, there is 
always a region in the left-half of the t-plane within which the expansion of C(t) in powers 
of u converges when — tt < < 0; but as approaches — 7r_ this region again shrinks to the 
point i on the cut. Recall that C(t) is the analytic continuation of C(t) from the right-half 
to the left-half of the t-plane. As long as |0| < tt/2 we can express C(t) as a convergent 
power series in u for all u inside the unit circle; in analogy with Eq. (|46|) we have 



<?(*) = XX(**K> ( 58 ) 

n=0 

where 

c n m) = (4>\C n (-t^H)\4,). (59) 

Recalling that = K\ip) and that K is the antiunitary time-reversal operator, so that 
K^K = 1 and ((b\K^)\a) = [(b\(K^\a))}* , we infer that for |0| < tt/2 we have 

CniU) = c n(U)- (60) 

This identity is hardly surprising since C(t) and C(t) represent the same (multivalued) 
function of t, and as long as |0| < tt/2 there is a common region in the lower-half of the 
t-plane within which the expansions of both C(t) and C(t) in powers of u converge. However, 
once crosses the line |0| = tt/2 into the region |0| > tt/2 the characters of Cnit^) and c n {tct>) 
differ since there is no common region in the t-plane within which the expansions of C(t) 
and C{t) both converge. On the cut along the negative imaginary axis in the t^-plane, c n (t n ) 
and c n (t- v ) are different branches of the same multivalued function of t^. 

An alternative form for C n (z), from which we may deduce the behavior of the coefficients 
c n{tij>) f° r large n, can be obtained using a standard expansion of the confluent hypergeometric 



function in terms of Bessel functions. []T6| We find that 



C n {z) = -2J2 A m (n)B rn+1 (z), n>l, (61) 



m=0 

where 



B m (z) = (-z/2n) m ' 2 J m (V^8^), (62) 
with J n (z) the regular Bessel function, and where A (n) = 1, Ai(n) = 0, A 2 (n) = 1, and 

A m+1 {n) = An-i(n) - [2n/(m + l)]A m _ 2 (n), m > 2. (63) 
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For large n and for z fixed, real and negative (so that \/—z is real and positive) we have, 
using the asymptotic form of J m (x) for large x, 

( _ z 2rn-l \l/4 j 1 

B ™{z) ~ ^ 22m+ i 7r2n2m+ i J cos(V=8^ - -mvr - -vr), n ~ oo. (64) 

It follows that for n ~ oo, and z fixed, real and negative, the first term in the series on the 
right side of Eq. ( |6"I| ) dominates, and we have 

/-2z\ 1/4 , 

C nW~ -ri cos(V-8n2 + 7r/4), n ~ oo. (65) 
V 7i n J 

Of course, —2 is not a fixed number, but rather is the operator t^H, which has a spectrum 
consisting of a discrete set of negative eigenvalues and a continuum of positive eigenvalues. 
From Eqs. ([12]) and fl47D we have 



Cra (t/>) — Cbd,n(^) + C CO nt,n (£</>), (66) 

where, using Eq. (|65|) , 

Cbd,n{t<f>) = X! IV'bd,m| 2 Cn(— t<f>Ebd,m) (67) 
m 

2 \ 1/4 / 

XI l^bd,m| 2 (^-Sbd,m) 1/4 cos(y / 8n^£' bdim + vr/4), (68) 



7r 2 n 3 



and 



Ccont,n(^) = / d'k |^(k)| 2 C„(-t^) (69) 

2 X 1 / 4 



rr 2 n 3 



J d 3 k (^S fc ) 1/4 | | 2 cos(^/8nt^ fc + 7r/4)e-" fc , (70) 



where in the last step we inserted an unobtrusive factor of e _,?A: (with 77 positive but in- 
finitesimal) to ensure convergence at a later stage. We choose the phase of E^m to be n, 
rather than — n, i.e. we write E^a,™, = e J7r | E^m \ , since the bound state poles are reached 
from the upper "physical" edge of the unitarity cut by following a path in the upper-half 
of the energy plane. Note that while J m (\/— 8nz) has a branch point singularity at z = 0, 
due to the square root in the argument, B m (z) does not have a branch point at z = since 
J m (x) is proportional to x m for x ~ 0. On the other hand, if we take the liberty of using the 
asymptotic form of B m (z) — see Eq. (0) — for complex values of z, we infer that B m (z) is 
not single-valued when <p is varied, with \z\ held fixed, around a closed loop from —71 to tv, 
and indeed we know that c n {t^) is not single-valued when <p is varied from — 71 to 71. 

Since -Ebd,™ < the bound-state terms on the right side of Eq. ( |68|) explode exponentially 
for n ~ 00, unless = ±n. Hence Cbd,n(^) also explodes exponentially, unless <fi = ±n 
(excluding a potential with a Coulomb tail, a case considered below.) This singular behavior 
is related to the fact that Cbd(£) has an essential singularity at t = 00, i.e. at u — 1. As t 
approaches 00 in the lower-half t-plane, Cbd(t) explodes exponentially, and unless <fi = ±7r 
there is a nonvanishing sector of the lower-half t-plane within which t approaches 00 and, 
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concomitantly, u approaches unity within the unit circle (this sector is the full half-plane if 
|0| = 0+). Therefore, unless = ±7r, the power series in u converges at u = 1 only if we 
exclude from \ip) the bound-state eigenvectors of H. Thus we replace \ip) by PcontlV') where 
-Pcont = 1 - -Pbd and where P con t and Pbd are projection operators with P bd defined as 

Pbd = X! IXbd,m) (Xbd,m|> (71) 
m 

where, if the potential is short-range (i.e. no Coulomb tail), the sum is over all bound states. 
The omission of these bound states does not affect the rate for a continuum process (if the 
potential is short-range), but it may affect the energy shift of the system. Note that while 
the compound states also give rise to an essential singularity, they need not be omitted; the 
sector of the t-plane in which they explode corresponds to letting u approach unity from 
without the unit circle, provided that |0| < H . 

To obtain the large-n behavior of c con t,n(^), we change variables from k to ^[n on 
the right side of Eq. (|70|) . Let us first treat the case of a short-range potential, and for 
simplicity let us assume that the angular momentum quantum number / is zero. We can 
factor \ip(k/ y/n)\ 2 out of the integral as |-0(O)| 2 , and we find that 

, , 15tt //i\ 3/2 |^(0)| 2 , x 

CcontAU) ~ 2972 [rj n ~°° (? ) 

where we used 

3/2 / \ 3 / 2 

/ d*k (t+Etf" cos( ^E k + 7r/4)e-^ = I . (73) 

We pause for a few remarks. First, while the integral on the left side of Eq. (^) is formally 
defined only for t^ real and positive, both sides of Eq. ( |72|) are analytic everywhere in the 
finite t^-plane cut along the negative real axis. Second, if I ^ we must modify Eq. (|72|) by 
including a factor proportional to 1/n 21 . Finally, note that the large-n behavior of c con t,n(^) 
is determined by small k, i.e. by the continuum eigenvalues of H close to threshold. Since 
the large-t behavior of C con t(*) is also determined by the continuum eigenvalues close to 
threshold, we infer that the large-t behavior of C con t(£) is determined by those terms with 
large n in the power series on the right side of Eq. (f46|). 

It follows that if the potential is short-range, c cont)n (/f:0) decreases as n~( 2l+5 ^ 2 as n in- 
creases. We now consider the case where the potential has an attractive Coulomb tail. 
We must first reconsider Cbd,n(^)- There are an infinite number of Rydberg bound states 
converging to threshold for which |8nt^i?bd,m| < 1 no matter how large is n. Such states 
do not yield divergent terms on the right side of Eq. (^3) and, as noted earlier, cannot be 
distinguished from continuum states converging to threshold; therefore, they should not be 
projected out of the wavepacket when ^ ±n. We write -Ebd.m = —Z 2 e 2 /[2(m*) 2 a ], where 
m* — m — 8, with 5 the quantum defect, which is roughly independent of m for m 3> 1. We 
must retain at least those bound states for which m* is roughly greater than or of the order 
of \fn. If m* > 1 we can relate l^bd,™! 2 to |-0(k)| 2 as follows. The wavepacket \ip) con- 
tains Rydberg bound states with population \iphd,m\ 2 dm* in the interval (m*,m* + dm*), and 
continuum states with population |?/>(k)| 2 ci 3 fc in the interval d 3 k centered at k. Just above 
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threshold, where kao <C Z, we can replace |?/>(k)| 2 by 27r7|-?/>(0)| 2 , where we used Eq. (|IDD 
noting that e~ 2nl <C 1. Since dE hd ^ m / dm* = Z 2 e 2 /[(m*) 3 a ] and dE k /dk = k/fjL, and since 
the differential populations with respect to energy must be the same for bound states just 
below, and continuum states just above, threshold, we have, writing d 3 k = 4irk 2 dk, 

\A d , m \ 2 =(^^)mO)\ 2 . (74) 



It follows from Eqs. flBlf ) and (|74]) that, with on the upper edge of the cut along the 
negative real axis of the t^-plane, and with E hd ^ m = e l7T \E hd ^ m \, 

/ 2 \V4 oo (f jp U/4 



c hd , n (t<t>) ~ 8n 2 Z*\m\ 2 {^) E ' (J;^ cos^Snt^m + vr/4), (75) 

where we have extended the lower limit of the sum to 1 since the leading n-behavior of 
the sum is relatively insensitive to the lower limit. In Appendix D we analyse the large- n 
behavior Cbd,n(^) and c con t,n(t/>) for all complex t§. To summarise the results of Appendix 
D, when the potential has an attractive Coulomb tail we find that 

i f 3ir 2 Z^(0)\ 2 \ 

Cn{U) ~ — , -7T < < 71, (76) 



n 2 \ 2 5 /2 ao t 



for tj. on the first sheet, and 



'37r 2 Z/i|^(0)| 



Cn(^) 2 ) ' 7T < < 37T, (77) 



n 2 



2 5 /2 ao t 



for ttj, on the second sheet. Hence, c n (t<p) has not only a branch cut on the negative real 
t^-axis, but also a discontinuity which amounts to a sign reversal. 

4 Green function 

The inclusive rate at which a continuous stationary or quasistationary process occurs if 
£ is the positive energy of the system is — 21m R(£ + irj) where i] is positive but infinitesimal 
and where R(E) is a Green function matrix element of the form 

R(E) = (^\G(EM. (78) 

The real part of R(S + irj) includes the energy shift of the system. For example, if the 
system consists of an infinitely heavy particle (at rest) and a light particle that is incident 
from infinity in the unperturbed state represented by |\?o)> we obtain the rate for the light 
particle to scatter by putting = H^l^o) where W is the interparticle potential. As another 
example, consider a system that consists of an atom, initially bound in the unperturbed state 
represented by \^q); if this atom is exposed to weak monochromatic radiation, we obtain the 
rate for the atom to decay by putting \ip) = V + \^q) where V + is the one-photon absorption 
operator. 
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Provided that E lies in the upper half of the complex i?-plane we can represent R(E) by 
the integral 



POO 

R(E) = -i dt e lEt C{t) 



Jo 



(79) 



roc 

= -it / dr e t{t » E)T C{t) 



Jo 



(80) 



Since C(t) is bounded throughout the sector < arg (r) < S we can rotate the contour 
of integration into this sector and extend the integral representation to a sector in the 
lower-right quadrant of the .E-plane. To analytically continue R(E) into a sector of the 
_E-plane that includes the negative real E-axis we first recall that both C TCS (t) and Cb g (t) 
are bounded throughout the lower right quadrant of the r-plane, but C^(t) is unbounded in 
this quadrant. Therefore, provided that we project out the bound states, and choose E to lie 
temporarily on the left edge of the positive imaginary S-axis, we can rotate the integration 
contour into the lower-right quadrant of the r-plane so that it runs along the right edge of the 
negative imaginary r-axis. Subsequently, we can move E from the upper- to the lower-left 
quadrant of the E-plane. Although we subtracted the bound-state contribution to R(E), 
this contribution is analytic, having only poles on the negative real E-axis. Hence, we have 
analytically continued R(E) from the upper- to the lower-left quadrant of the E-p\ane. 

To deal with all values of E in the lower half of the complex E-pl&ne we introduce 
R(E) = (ifj\G(E)\ifj) , where we recall that = K\ip). For Im E < we can represent 
R(E) by the integral 



Since C(t) is bounded throughout the sector 7r — S < arg (r) < ir we can rotate the contour 
of integration into this sector, and again extend the validity of the integral representation, 
now to a sector in the upper-right quadrant of the .E-plane. To analytically continue R(E) 
into the upper- left quadrant of the i?-plane we proceed as for R(E); we project out the bound 
states, choose E to lie temporarily on the left edge of the negative imaginary i?-axis, rotate 
the integration contour through the lower left quadrant of the r-plane, so that it runs along 
the left edge of the negative imaginary r-axis, and subsequently move E into the upper-left 
quadrant. 

Since C(t) and C(t) are continuous in the lower-half t-plane, the integral representations 
of R(E) and R(E) are identical when both contours are chosen to be the negative imaginary 
r-axis. Hence R(E) and R(E) are identical in the lower-left quadrant of the -E-plane; it 
follows that R(E) is the analytic continuation of R(E) from the upper part of the i?-plane 
to the lower part, on the first sheet of the Riemann energy surface cut from to oo along 
a line parallel to the integration contour of either R(E) or R(E). At the cut, R(E) and 
R(E) are different branches of the same function. As noted above, the integration contours 
of R(E) and R(E), respectively, can be rotated into sectors of the lower- and upper-right 
quadrants of the E-pl≠ hence the cut can be rotated into these sectors, thereby extending 
the region of the Riemann energy surface covered by the integral representations. However, 
these integral representations do not permit the entire Riemann surface to be covered. As we 
see below, the entire surface can be covered by the series representation. The wider region 




(81) 
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of analyticity afforded by the series representation is a reward for performing the integration 
over time. 

Using the expansion of C(t) in powers of u we have 

R(E) = -it / dr £ c n m)u n e^ T . (82) 

Rotating the integration contour to the line arg (r) = <p, interchanging sum and integral, 
and recalling that 



00 



l n (a) = / dr e iaT (— , (83) 

JO \T — lJ 

we arrive at the series representation 

00 

R{E) = -it+ c n m)l n (Et^), (84) 

n=0 

where for the moment \<ft\ < it/2. We obtain the same series representation of R(E) after 
expanding C(t) in powers of u, and rotating the integration contour to the line arg (r) = 7i+<f). 
The resolvent itself has the series representation 

PcontG(E) = -it P cont / dr e l{Eto)T U(t) (85) 

Jo 

00 

-it^Y. IniEt^P^Cni-t+H), (86) 

n=0 

as discussed in the Introduction. Since R(E) is defined by an integral whose contour is the 
positive real r-axis, this contour must be mapped into the unit circle, as it is if <fi lies in the 
range < <fi < 7r/2; hence the series representation of R(E) is defined for (f) in this range. 
On the other hand, R(E) is defined by an integral whose contour is the negative real r-axis, 
which is mapped into the unit circle if <fi lies in the range —71/ 2 < <fi < 0; hence the series 
representation of R(E) is defined for </> in this range. It follows that when E is real and 
positive, the physical branch is specified by < <fi < 7r/2, whereas the unphysical branch is 
specified by — 7r/2 < <p < 0. 

We now prove the useful result that we can pass (discontinuously) from one branch of 
R(E) to the other branch, i.e. to R(E), by simply changing the sign of the phase <fi. We 
first observe that, according to Eq. ([43]), changing the sign of <fi and complex conjugating 
u amounts to changing r into — r*, or, equivalently, t into —t*. It follows from Eqs. fl4"S| ) 
and ([5j]) that complex conjugating C(t) and changing the sign of yields C(—t*), which, 
using Eq. (|3"6"D, is the same as [£?(£)]*. Hence, simply changing the sign of 4>, without 
complex conjugation, turns C(t) into C(t). Furthermore, changing the sign of <j), from 
positive to negative, say, transforms the contour of integration for R(E) discontinuously into 
the contour for R(E); if we vary <p continuously from a positive to a negative value the 
integration contour — the broken line in Fig. 3 — starts as a path in the upper semicircle, it 
moves upwards towards to the boundary of this semicircle, and upon reaching this boundary 
jumps discontinuously to the boundary of the lower semicircle, and from there to a path 
inside the lower semicircle. Thus we have arrived at the result we set out to prove. 
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4> = 7t 




Figure 4: The physical sheet of the Riemann energy surface cut along the negative real axis. 
Bound state poles (x) lie on the upper edge of the cut and resonance poles (*) lie in the 
lower-half plane. The series representation of the resolvent applies on this sheet if we choose 
the phase parameter to be = 7r. 
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An important task that remains is the analytic continuation of the series representation 
of R(E) from |0| < tt/2 to \(f)\ > tt/2. To carry out this task we must analyse the integrals 
X n (a), defined by Eq. Q8lf). These integrals are not formally defined when Im a < since 
the integrands do not vanish for r ~ oo. We can analytically continue T n (a) to the region 
Im a < by rotating the integration contour into the upper-half r-plane. However, if n > 1 
the integrand is singular along the positive imaginary r-axis, and therefore the (positive) 
angle through which we rotate the contour must be less than tt/2, which only permits analytic 
continuation of Z n (a) to the lower-right quadrant of the complex a-plane. To analytically 
X n (a) continue to the lower-left quadrant of the a-plane we move a temporarily to the upper 
edge of the negative real axis in the a-plane, rotate the integration contour through the angle 
— 7r, so that it runs along the lower edge of the negative imaginary r-axis, and subsequently 
move a into the lower-half of the a-plane. However, when n > 1, X n (a) has a branch 
cut extending from to oo in the a-plane since we cannot move a continuously around 
a closed loop by continuously rotating the contour of the integral representation of X n (a). 
The direction of this cut is a matter of convention, and that convention has already been 
established. Inspection of Eq. ([86]) reveals that the choice of the cut for X n (a) determines the 
cut in the Emplane for the series representation of R(E), and we want that to be the same as 
for the integral representation of R(E). Thus we draw the cut for X n (a) along the positive 
real axis in the a-plane, which implies, according to Eq. (^), that the series representation of 
R(E) has a cut in the E-plane along the line arg E = —(f). This conforms to our specification 
that R(E) and R(E), respectively, are analytic in the upper and lower halves of the E-plane, 
and have integral representations defined in the ranges < < tt/2 and —tt/2 < < 0. 
We now see how the series can represent R(E) over the entire Riemann energy surface. If 
we let increase from to tt_, the branch cut swings from the positive real E-axis, down 
past the resonance poles in the lower-half E-plane, to the lower edge of the negative real 
E-axis. Hence, fixing to be tt_, we see that R(E) is defined over an entire energy plane, 
and its poles are the bound-state poles on the upper edge of the cut and the resonance poles 
in the lower half of the energy-plane. See Fig. 4. This energy plane is the physical sheet 
of the Riemann surface, since R{E) is the branch corresponding to outgoing-wave boundary 
conditions when E is real and positive. Similarly, if we let derease from to — 7T_, the 
branch cut swings up past the resonance poles in the upper-half .E-plane, to the upper edge 
of the negative real E-axis. Hence, fixing to be — vr_, we see that R{E) is also defined over 
an entire energy plane, and it has bound-state poles on the lower edge of the cut; but its 
resonance poles are conjugate to those of R(E), as shown in Fig. 5, and the energy plane 
is the unphysical sheet of the Riemann surface since R(E) is the branch corresponding to 
ingoing-wave boundary conditions when E is real and positive. 

To determine the rate of convergence of the series representation of R(E) we need to 
analyse the behavior of X n (a) for large n. The integrand on the right side of Eq. (|83| ) 
oscillates rapidly when either n or a is large. To treat large n it is useful to change to the 
variable 9 = tan _1 (r). We have 



We start by assuming that a is real and positive and that n > a/2. There is a point of 




(87) 
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Figure 5: The unphysical sheet of the Riemann energy surface cut along the negative real 
axis. Bound state poles (x) lie on the lower edge of the cut and resonance poles (*) lie in 
the upper-half plane. The series representation of the resolvent applies on this sheet if we 
choose the phase parameter to be = — n. 
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stationary phase at 9 = 6q, where 



cos 2 # = a/(2n), (88) 
and if n » a/2 we have, using the method of stationary phase, 

/ 2 \ 1/2 

X n {a) ~ (-1)" sec 2 9 e ins(e )+i^ (g9) 

\n\s (0 )\J 

where s" (9) is the second derivative with respect to 9 of 

s(6) = {a /n) tan 9 - 29. (90) 

Noting that s"(9 ) = 4tan6 l and that 

tan^o « (2n/a) 1/2 (l - a/ An) 

we find that ^ 

X n ( a ) ~ f ^) e ^+^ 4 , a > 0. (91) 



a" 



Evidently, when a is real and positive X n (a) grows as n 1//4 with increasing n. When a is real 
and negative, Re X n (a) = 0; this follows upon rotating the contour of integration on the 
right-side of Eq. (p3|) from the real positive axis to the negative imaginary axis. For complex 
values of a, in the plane cut along the positive real axis, the point of stationary phase 
becomes a saddle point, and the expression on the right-side of Eq. fl9T|), when analytically 
continued, becomes exponentially small, i.e. it decreases exponentially with increasing n 
as e~^ m v^<™. However, we must include not just the contribution from the saddle point, 
but also the "surface" contribution, at the lower limit of the integral X n (a); integration by 
parts yields the following asymptotic expansion in powers of 1/n, applicable as long as the 
saddle-point contribution is negligible: 

n J^o n 

where fj(a) is a polynomial of degree j with real coefficients. For example, fo(a) = 1/2, 
/ x (a) = a/4, and / 2 (o) = (a 2 - 2)/8. 

Suppose that E is real and positive. Let us set = 0± in the series representation of 
R(E). The large-n behavior of the factor X n {Et^) is given by Eq. fl91|) . These factors, while 
oscillating as n increases, also weakly magnify each term of the series for R(E). However, 
since the c n {t^) decrease more rapidly than theX n (£^0) increase, the series converges, albeit 
logarithmically. We are now at the stage where we can work through some examples which 
support the validity of the series representations of G(E) and R(E), i.e. Eqs. (|86D and ([84]). 
We do this in Appendix B. 

If is chosen to be in the range < \<f>\ < tt, the saddle-point contribution to X n {t^,E) 
can be neglected for n sufficiently large, i.e. for \/8ntoE sin(|0/2|) ^> 1. In this case we can 
use the asymptotic expansion of Eq. ([92]) to express, in a rather useful form, the contribution 



28 



of the higher terms to the series representation of G(E). Combining Eqs. flSBD, (|57J|), and 
(|54l) gives 

PoontG(E) = P cont e-^ H Q + tt^H)S(E, 0)) , (93) 
where S^-E, 0) is the infinite sum 

oo 

0) = E 2„(VE) iFi(l - n, 2, -2fy#)- ( 94 ) 

n=l 

Using Eqs. fl9"2] ) and (Q), and neglecting the contribution from the saddle point, gives, after 
interchanging the order of the sums, 

oo oo / i\n 

S(E,<f>)~-iY, /i-i(^)E ^^(1-^2,-2^). (95) 

j=l n=l nJ 

We now use the Watson transform: 

00 1 r F(n) 

Y (-l) n F(n) = - dn p^-, (96) 

^ v ' v ' 2iJv sin(nvr)' v ; 

where T is a counterclockwise contour which encloses the section of the real positive n-axis 
extending from a point in the interval (0,1) to infinity. We distort V into a line, T', parallel 
to the imaginary n-axis, with < Re n < 1, a region where it is permissible to use the 
integral representation 

sin(n7r) r 1 , (\ — u\ n . . 

Ml-n^z) = — i — -/ due zu { . 97 

nn Jo V u J 

It follows that 



i=i 



■it 



Define uo = ln[(l — «)/«]. If w < 1/2 we have u > 0, while if w > 1/2 we have < 0. 
Now, if M > 1/2 we can bend the contour V back into V; we obtain zero since T no longer 
encloses any poles. If u < 1/2 we can bend V into a contour T" that encloses the whole 
of the negative real n-axis and, in addition, a section < n < n < 1 of the positive real 
n-axis. We use, integrating by parts in the first step, 

dn — — = — dn 99 

n 3+v ]\ Jt" n 

= U -( dn e — + J dn 6 —), (100) 

J I \J — oo ii Jn^—iT] IL / 



where 7/ is positive but infinitesimal. Transforming the integration variable from n to —n/cu, 
and noting that the exponential integral, E\(z), has a logarithmic branch point at z — 0, 
with a branch cut customarily drawn along the negative real z-axis, we obtain 

\ dn = Ei(—ujn + irj) — Ei(—ujn — irj) (101) 

Jt" n 

= -2m. (102) 
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Putting all this together, we obtain 



S(£,0) ~zfj tt^h3 £ /2 due- 2t * Hu 



In 



■it 



it, 



(103) 



The evolution operator e~ 2t 4> Hu acts for at most one unit of time, and can be calculated 
rapidly (using a Pade approximant — see Appendix C). The first few terms — say the first 
iV — on the right side of Eq. (OT) should be evaluated exactly, and accordingly the first iV 
terms built from the right side of Eq. ( |92|) should be subtracted from S(E, 0), with N chosen 
so that for n > N the following conditions are satisfied: (i) The asymptotic expansion of 
Eq. fl92|) converges when only a small number of terms are included, and (ii) the saddle-point 
contribution to I n {t^E) can be neglected. 



5 Conclusion 



We have derived a series representation of the resolvent G(E) which incorporates a com- 
plex unit of time, toe 1 ^, whose phase specifies the branch. This permits a real basis to 
be used to construct a matrix representation of the Hamiltonian. The foundation of our 
approach is the analytic continuation of the temporal correlation function in the complex 
time-plane. Using analytic continuation we can extrapolate the correlation function from a 
relatively small interval, of duration to, to an interval of asymptotically large duration. Since 
much of the dynamics take place during a time that, typically, is comparable to to, it should 
be possible to employ a basis efficiently; it is unnecessary to span spatial regions that are 
reached by the particle(s) at times very much larger than t Q , where the wavefunction takes 
on its asymptotic form. 
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Appendix A: Equivalence of Two Forms of the Integral T n {a) 

In this appendix we show that Eqs. (|) and @ are equivalent definitions of T n (a). To do 
this we use a standard contour integral representation of the Laguerre polynomial: 

L ( nUx) = ^- I e~ x <dC 
2m J 

where the contour runs counterclockwise and encloses the point ( = 0. Substituting this 
representation into the right side of Eq. (H) gives 

TJ Ett) = -\ f « (i±i)" f M e-™«0 - L^) . (105) 

The integrals of Eq. (|105|) are defined as long as (i) E is not real and positive and (ii) 
Re t(j,{l + 2() > 0. This second condition implies that the closed contour excludes the point 



i + C 



(104) 
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( = —1/2. The first condition can be dropped after integrating over H. Performing the 
integration over H gives, with a = Et^, 

*■<«> = s/« (^)"(^ (I+2C ^i- 0(1 + 2C )i + ^iVc))- (106) 

where the branch cut of Ei(z), the exponential integral, is along the negative real z-axis; the 



contour of integration in Eq. ( |106|) excludes the branch point at £ = —1/2 and is drawn to 
exclude the cut along the line arg a(l + 2Q = 0. Changing variables from ( to r = 1 + 2( 
we obtain 

IJa) = yJ' 1t (Bt)" » + wtt) • < 107) 

where the new (still counterclockwise) contour excludes both the point r = and the branch 
cut along the line arg ar = 0, but includes the point r = 1. We now deform this contour so 
that it wraps around the branch cut and runs along the circle at infinity. Since e z Ei(z) ~ 
(1/z) — (1/z 2 ) for z ~ oo the integral over the circle at infinity vanishes. The integral around 
the branch cut can be simplified after noting that as z moves from the upper to the lower 
edge of the negative real z-axis, Ei(z) increases by 2ni. It follows that 



-oo 



I n (a) = i I dr e- ar ( ) . (108) 



o 



By rotating the integration contour through 90°, from the negative real r-axis to to the 
negative imaginary r-axis, we reproduce Eq. (|j). 

Appendix B: Analytical Tests 

In this appendix we give two examples where we use the series representations given 
by Eqs. (^) and (|84]) to establish some known formal results. In our first example we 
use Eq. fl36|) to derive the standard expression 1/(E — H) for G(E) when the potential is 
arbitrary; but we assume that the part of G(E) which accounts for propagation far off the 
energy shell is relatively unimportant. We choose both t an d <p to be small, such that 
Et is significantly smaller than 1 and \<f>\ <C n. Hence the exponential e l V 8t 'i> En is almost 
undamped and is slowly oscillating. Therefore a large number of terms contribute to the 
series on the right side of Eq. (^) and the sum is determined primarily by large-n terms. 
Consequently, we can substitute the asymptotic forms for C n (— t^H) and X n {Et^) — see 
Eqs. ( |65|) and (^l|) — into the series of Eq. (|^) to give 



G(£)^-^(-j £ (_) e ^W4 cos( ^ + ff/4) . (109) 



The cosine on the right side of Eq. ( |109| ) has two parts, one co-oscillating and the other 
counter-oscillating relative to the pre-exponential. Neglecting the co-oscillating part (its 
contribution averages to zero) gives 
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Replacing the sum over n by an integral over n, extending the lower limit of integration to 
zero, and changing variables to y = \/8En, we obtain 



G(E) « -i^ilT Cdy^^ (111) 

(112) 



2E \E, 

1 fH\^ 4 1 



2E \EJ 



H/E) 



where we used J °° dy e ty = i. Recall that when n is large, C n (— t^H) contributes to the 
correlation amplitude C(t) at large times. Therefore terms with large n in the series repre- 
sentation of G(E) contribute primarily to the near-energy-conserving part of G(E); terms 
with small n, which we have neglected, contribute to the non-energy-conserving part of G(E). 
Hence it is consistent for us to replace the overall factor (H/E) 1 / 4 by unity, and the term 



(1 - y/H/E) by (1 - H/E)/2; this yields G(E) = l/(E - H). 

Now we discuss our second example; we derive from Eq. fl84|) the threshold behavior of 
the rate — 21m R(E) for E ~ 0. Again we need focus only on large- n terms on the right side 
of Eq. (0). For simplicity we assume that the angular momentum quantum number / is 
zero. Assuming also that the potential is short-range (i.e. no Coulomb tail) it follows from 
Eqs. (||D, (]72|), and (0) that for E ~ 



iKf ..\3/2 \ oo 

R(E) ~ -^(0)| 2 1/9 lb ™ y ^ . (113) 

^^(2^3)1/4^^ „9/4 



Replacing the sum over n by an integral over n with lower limit no ^> 1, and changing 
variables to y where y 2 = y8En, we obtain 

r(e) ~ -iV^iv(o)i 2 (7r M ) 3 / 2 (^^nj y w -j^g-, (H4) 

where y = (SEno) 1 ^. As E approaches zero, so does yo, and the integrand becomes highly 
singular at the lower limit. However, before taking the limit E — > we can regularize the 
integral by integrating by parts three times, and by discarding the surface terms at the 
upper and lower limits; the surface terms are exponentially small at the upper limit and are 
of order l/fyot^ 4 ) or less at the lower limit, and therefore negligible if no is sufficiently large 

that 8En t » 1. Note that the exponential e ty of Eq. ( p.!4|) restricts y to values less 
than or of order <p~ l / 2 tQ 1//4 and hence l/(yo£c/ 4 ) cannot be less than a number of order 1 / 2 ; 
therefore we require <p ir. After integrating by parts three times we can let the lower limit 
of the integral be zero, and we arrive at 

R{E) ~ -^4 /4 |^(0)| 2 (7r/i) 3 / 2 2 7 / 2 e i7r / 4 J~ 'dye^V**. (115) 

Performing the integration over y yields, for E ~ 0, the threshold law 

- 21m R(E) ~ 2 7 / 2 |^(0)| 2 ttV/ 2 v / ^, (116) 
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an expression which of course is independent of t^. This is the correct threshold law if I = 
and if the potential is short-range. Indeed, as a check we can use Eqs. ( [TED and (|79|) to write, 
for E ~ 0, 

poo piEt 

- 2Im ME) ~ 2 5 / 2 (7r/x) 3/2 |^(0)| 2 Re e"^ 4 / dt . (117) 

jo (t + 7/)^ 

where we have introduced a positive but infinitesimal term r\ in the denominator to regularize 
the integrand at t = 0. Integrating the right side of Eq. (|117|) by parts (once), and evaluating 
the resulting integral, we reproduce exactly the right side of Eq. ( |116| ). For I ^ we must 
modify Eq. ( |116| ) by including a factor proportional to E on the right side. 

If the potential has a Coulomb tail it follows from Eqs. fl53|), ([75]), and ([JT]) that 

fl(s > - w°>l 2 (^) £ (118) 

Replacing the sum over n by an integral over y, where again y 2 = \^8En, yields 

R(E) ~ 12^(0)1^ M/I / d y -^i ( 119 ) 



and, after regularizing the integral by integrating by parts twice, and subsequently integrat- 
ing over y, we obtain 

R(E)~-i8ir 3 (Zfi/a )\iP(0)\ 2 , (120) 
which gives the energy-independent threshold law 

- 21m R(E) ~ 16vr 3 (Z/i/a ) |V>(0) | 2 . (121) 
As a check we can use Eqs. ([H]) and (|79|) to write, for E ~ 0, 

-2Imi2(E) ~ 327r 2 (Z/i/a ) |^(0) | 2 Im / dt ; (122) 

Jo t 

since dt sm(Et)/t = tt/2 we reproduce Eq. ( p.21| ). 
Appendix C: Algorithms 

In this appendix we describe some algorithms for the numerical implementation of the 
series representations of G(E) and R(E). We consider the evaluation of (i) C n (z), with z a 
linear operator, (ii) X n (a) for all complex numbers a, and (iii) e z , with z a linear operator. 

(i) The Operator C n (z) 
Recall that C n (z) = Q n (z)e z , where Q n {z) is a polynomial of degree n in z: 

Q n {z) = {2z) 1 F 1 {l-n,2,-2z), n > 1, (123) 

with Qo(z) = 1. Using a standard recurrence relation for the confluent hypergeometric 
function ]n| we obtain the recurrence relation 

(n + l)Q n+1 (z)=2(n + z)Q n (z)-(n-l)Q n „ 1 (z), n>l, (124) 
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which can be started using Qi(z) = 2z. Alternatively, we can formulate a backward recur- 
rence relation for the function B m (z) of Eq. (|62"D using a standard recurrence relation for the 
Bessel function. [TB| We find that 

B m -i(z) = --[mB m (z) - 2nB m+1 (z)], m > 1. (125) 

Since B m (z) decreases exponentially as m increases, this backward recurrence relation can 
be started by putting B m+ i(z) = and B m {z) = 1 for a large value of m, and by using the 
identity |16| 

00 /— 2n\ m 

B (z)+2J2 B 2m (z) = l (126) 

m=l ^ Z ' 

to correctly renormalize the functions. After computing the B m (z) we can use Eq. (pH) to 
calculate C n (z). 

(ii) The Integrals X n (a) 



To derive a recurrence relation for the integrals X n (a), with Im a > 0, we start from Eq. 
Using 

cos 2 9 = (2 + e 2i9 + e- 2ie )/4, (127) 

we see that 

21 n {a) -l n -i{a) -X n+1 (a) = 4(-l) n d6 e iat ™ e ~ 2me . (128) 

Jo 

Integrating by parts on the right side of Eq. (|128|) , and discarding the surface term at the 
upper limit 9 = tt/2 (it vanishes if Im a > 0) yields 

21 n {a)-l n ^{a)-l n+1 (a) = (-l) n+1 (^-J + (-1)™ (^-J J d6 sec 2 6e iat ™ 8 - 2m0 . (129) 

Recognizing that the integral on the right side of Eq. ( |129| ) is proportional to T n (a) we arrive 
at the recursion formula 

ln + M = (-l) n (-) +2(l--) X n {a) - J n _i(a), n > 1, (130) 

which can be started using 



n J \ n 



X (a) = i/a, (131) 
2i(a) = X (a) + 2ie~ a E!(-a), (132) 

where -E , 1 (^) is the exponential integral. As long as a is positive (with an infinitesimal positive 
imaginary part) X n (a) increases (albeit weakly) as n increases, and this forward recurrence 
relation is stable. However, for large a we have (see below) X n (a) ~ (— l) n (i/a) and therefore 
there is cancellation between the first and second terms on the right side of Eq. (130). Thus 
we consider separately, below, the case a^> n. 

Note that E\(z) has a logarithmic branch point singularity at z — 0, and this singularity, 
which first appears in Eq. ( |132| ) , is propagated by the recursion formula of Eq. ( |130| ) so that all 
terms but the first on the right-side of Eq. ( |34"D have logarthmic branch points at E = [the 
branch cuts lie along the line arg (E) = —0]. In contrast, the exact R{E), which is a function 
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of the dimensionless variable Eto, has a square- root branch point at E = (if the potential 
is short-range). Presumably the converged sum on the right-side of Eq. (jS4j ) does have the 
correct branch-point behavior. In this regard it is helpful consider an example, namely the 
function f(z) = Y^Lq [ln(z) / 2} n / n\ . Any approximation to f(z) obtained by truncating the 
sum has a logarithmic branch point; but the infinite sum is f(z) = exp[ln(z)/2] = z 1 ^ 2 , 
which has no logarithmic branch point, only a square root branch point. 

If H »fiwe can represent X n (a) by an asymptotic power series in 1/a. Integrating by 
parts on the right side of Eq. (1531) yields 



- ^ M (0), (133) 

m=0 V0,/ 

where l^ m ^(r) is the m-th derivative with respect to r of 

Y^(r) = fI±lV '. (134) 



We have Y^(0) = (—1)™ and the derivatives can be calculated using the following relation 
(which can be obtained after some straightforward algebra) 

[m/2] , 

Yt +1 \r) = 2n^ ("^r^^^^O), (135) 

where [m/2] is m/2 rounded to the lowest integer. It follows from Eq. ( 135 ) that Y r [ m \0) 
is a real number multiplied by i m , and hence if a is real each term of the asymptotic series 
(|133|) is purely imaginary. In fact, the real part of 1 n {a) is exponentially small; to see this, 
rotate the contour of integration on the right side of Eq. ( |55| ) by 180°, counterclockwise. 
This yields X n (a)*, and since there is a pole at r = i it follows from Cauchy's theorem that 
the real part of Z n (a) is ni multiplied by the residue of this pole. Thereby we arrive at 

ReJ„(a) = e-(-)x; ^t^Tl (136) 

\aJ m\[n — 1 — m)\\n — my. 

= - 2ne £ (* + !)!(„ -*-!)! *! (137) 

= -27re~ a Lj l 1 2 1 (2a), Im a = 0, a>0 (138) 

Note again that if a is real and negative, Re T n {a) = 0. 

(iii) Exponentiation of a Linear Operator 
To evaluate e z we use the Pade approximant 

1 _|_ £ _L *L 
2 12 

which matches the power series expansion of e z through the term in z 4 . Since z is an 
operator it is convenient to factorize this expression so that we do not incur the additional 
computation of z 2 ; thus we write 

^ ^ (l-z/ Zl )(l-z/z 2 ) 



{i + z/z l ){i + z/ Z2 y 



(140) 



35 



where the roots z\ and zi are 

z 1:2 = -3 ± iV3. (141) 

If the Hamiltonian H is represented by the matrix H_ , constructed from a basis with overlap 
matrix B_, the time-evolution operator e~ %Ht becomes e~^— — )*, and we have (letting £ — ► 
-it) 

e ""' " R - (l M H l - + {tMm B - (k)fl '- + M^]. (142) 
This is an extension of the standard Cayley form of the time-evolution operator from third- 
to fifth-order; we have found the stability and accuracy of exponentiation to be significantly 
improved by this increase in order. While it is unnecessary to calculate 5 _1 here, note that 



since B_ is real, symmetric, and positive definite it has a Cholesky decomposition [I7j and 
hence its inverse can be calculated rapidly. 

Appendix D: Large- n form of c n (t^) when Coulomb Tail is Present 

To establish Eqs. fl76|) and Eqs. (|77|) for the large-n form of c n {t<p) when a Coulomb tail 
is present, we first show how to analytically continue Cbd,n (£</>) along a path in the upper-half 
£^-plane, from the upper edge of the cut to the positive real axis on the first sheet of the 
Riemann ^-surface. To this end we seek to replace the sum over m on the right side of 
Eq. (f74[) by an integral over m. However, the direct replacement of sum by integral is not 
entirely justified since the summand varies rapidly when n ~ oo, and the m-th and (m + l)-th 
terms of the sum may differ significantly. Nevertheless, we can accomplish the passage from 
sum to integral by making a small modification to the integrand. We first break the cosine 
in the summand into a sum of two exponentials: 

c hd , n (U) ~ 8tt 2 Z 3 |^(0)| 2 (-j^) [S+(U) + SL(£^)], (143) 

where 

2(m*a ) 3 ' U44j 



and we introduce the integrals 



(^-^bd.m) 1 ^ 4 i^/8nUE hd , m +i-K/4 



and 



I+iU) = / dm (M5) 
Jc+ 2(m*a y i 

UU)= dm [H hd ' m > — s , (146) 

v ^' Jc. 2(m*a ) 3 (l-e 2Mrm ) V ; 

where the contours C± run from m = 1 to oo along the upper edge of the positive real 
m-axis (and where m* is now the continuous variable m — S). We have introduced the factor 
(1 — e 2l7Tm ) into the denominator of the integrand of J_(£^), a factor which vanishes when m 
is an integer. The integral I+{t<f>) is formally defined for all t<f, on the first sheet, except on the 
cut, of the Riemann ^-surface, i.e. for — n < < it, since iJSntsE^ m = — w8n£j£' bdm | 
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and the exponential in the integrand of I+(t<f,) decays as m increases. Furthermore, due to 
this exponential, only the region in which | Snt^i^d,™ I is less than or of order unity contributes 
significantly to in this region m*, and hence m, are greater than or of the order of 

y/n, and so the integrand varies slowly with m. Consequently, when t$ is on the upper 
edge of the cut the integral and the sum S+(tj) are the same. The integral I-(t^) is 

formally defined for all on the second sheet reached by crossing the branch cut, i.e. for 
7i < (f) < 37T. For on the branch cut, the exponential in the numerator of the integrand 
of I-(tj,) oscillates, but it is undamped so we cannot directly replace the integral by the 
sum. However, we can express the integral as a sum by observing that the integrand of 
J_ (tc/,) has an infinite sequence of poles on the positive real m-axis, on the lower edge of the 
contour C_, at those points where m is a nonnegative integer. If we rotate C_ downwards, 
past the line of poles, we pick up the contributions from the poles according to Cauchy's 
residue theorem. The integral along the new, rotated, contour is negligible for the following 
reason: The exponential in the numerator of the integrand is now damped, and so the main 
contribution to the integral comes from the region where m is greater than or of the order 
of \fn\ but the exponential e 2lnm in the denominator is exponentially large since Im (m) is 
large and negative along the important segment of the rotated contour. It follows that for 
on the branch cut, and anywhere on the second sheet, I-{t^) is, when n ~ oo, just the 
sum of the residues of the poles multiplied by this is the same as S_{t^) when is on 

the cut. 

To evaluate I+(t^) we change variables from m to y — \J Zej (m*a^ 2 ). The contour 



of integration runs from y = to roughly y Ze/(a^ 2 ). Extending the upper limit to oo 
(which does not affect the leading n-dependence of the integral) and rotating the contour of 
integration downwards through an angle of 0/4, replacing y by ye - "^ 4 , gives 



-7r < (j) < TV. (147) 



\Z 2 a ) 

i ( 37T 1 / 2 /! N 

n 5/4 [2^Z 2 a t (f>/ 

A similar change of variables for 7_ (t^) gives, after extending the upper limit of integration 
to oo and rotating the contour of integration upwards through an angle of 0, 

= e i5e tti) d y y 4 5™—. (148) 

The integral I-(t$) is formally defined if (37r — 0)/4 > > (ir — 0)/4. We set equal to its 
minimum value (n + — <f>)/4. Let us, for the moment, choose t$ to be on the second sheet, so 
that > 7r; this implies that < 0. Due to the exponential in the numerator on the right 
side of Eq. ( |148j ), the main contribution to the integral comes from values of y less than or 
of order l/n 1//4 , and since < the exponential e 2l7r( - e 2t@Ze / a o / v 2 ) in the denominator on 
the right side of Eq. ( |148j ) is negligible. It follows that 
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[The right sides of Eqs. ( |147| ) and ( |149|) differ only by a sign.] Now let us choose t$ to be 
above the branch cut on the first sheet, so that < n; this implies G > 0. Again, the main 
contribution to the integral comes from values of y less than or of order l/n 1//4 , but since 
> the exponential e 2ln( - e 2lBZe / a o v 2 ) is now very large. Hence, as t$ crosses the cut from 
below to above, moving from the second sheet to the first, I-it$) vanishes through order 
1/n 5 / 4 . It follows that for t$ on the first sheet 

, , i /37r 2 Z/i|^(0)| 2 \ 1 , 1ro , 

c M ,„fa)~^( 2v^ f ; j- -*<*<'• (l50) 

but that Cbd,n{t(j>) vanishes through order 1/n 2 on the second sheet, i.e. when n < < 3n. 
For 7^ ±7r_ we must project out those bound states for which m* is roughly less than y/n. 

We now analyse c contjn when the potential has an attractive Coulomb tail. The 
integral over k on the left side of Eq. (|70|) is, strictly speaking, defined only when = 0, 
since cos(J8nt t j ) Ek + n/A) diverges as k increases unless t<f, is real and positive. However, as 
before we break the cosine into a sum of two exponentials. After performing the integration, 
we can analytically continue the separate integrals. Let us start with = to and move 
along a path in the upper-half plane to the point toe l7T ~ . We first change variables from k to 
y = (k/ yZ/i) 1 / 2 , which removes from the integrand the branch point arising from the factor 



A; 1 / 2 . Note that i^J&nt^Ek is equal to 2ie % ^/ 2 y 2 ^/nt^, and has a negative real part over the 
entire path. It follows that, after breaking the cosine into a sum of two exponentials, the 



integral over the term in exp (i J Snt^Ek + in/A) is well-defined over the entire path since its 
integrand decreases exponentially as y 2 increases. On the other hand, the integral over the 



term in exp^i^JSnt^Ek — in/A) is not well-defined; its integrand increases exponentially as 
y 2 increases. To analytically continue this second integral we proceed as before, and rotate 
the contour of y- integration downwards through an angle G where > 0/4 so that the 
exponential does not explode with increasing y 2 as we move t^. The first integral is 



\j d*k{UE k )^\m\ 2 e l ^^ +m/A = ^(^y 



1/4 







where 7 = Z/ (ao^^y 2 ) and where dVL is an element of solid angle containing k. We now 



change variables in the first integral from y to x = yn 1 ^ so that k = x 2 ^Jfi/n and 7 = 
(Z/aoX 2 )J n/fj,; for n ~ 00 we can neglect e _27r7 compared to 1, and we can replace |?/>(k)| 2 



by |-0(O)| 2 , so the first integral becomes 

Mt\mi(f^) ut i da r dxA **>yu = i(^ymf) (152) 

a n 5 / 4 \ 2 I J Jo n 5 / 4 I 2 n /4 ao t I v > 
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The second integral is, after changing variables from y to ye 



ie 

1 



x f^^e- 2 ^-- 2 "^, (153) 

where now k = [i 1 / 2 y 2 e~ 2ie and 7 = Ze 2lQ / (ao^^y 2 ) . We fix B to have its minimum value, 
i.e. B = 0/4+^/2, with r\ positive but infinitesimal and change variables from y to x = yn 1 ^. 
For n ~ 00 we have 

87i 2 Zfie~ i5e \tP(0)\ 2 /e"%\ 1/4 



- / d 3 k (t^^WlVV 8 ™^--/ 4 



/■oo g-2(i+?7)a; 2 v / io 

x / dxx 4 - (154) 







1 - e- 2 ^ 



where k = x 2 e~ 2 * e 1 //i/n and 7 = e 2%e [Z ' j a x 2 ) J nj /i. Note that the integrand has an infinite 
number of poles, accumulating at y = on the line arg (y) = —it/A, at those points where 
e -27T7 _ ^ ^ g j on g ag _ 7r < < 7]- we have — 7r/2 < 2B < 7r/2, so that the contour of 
integration lies apart from the line of poles, and, furthermore, Re 7 ^> 1 so that the term in 
e -2?r7 Qn ^ e right side of Eq. (|154j) is negligible; hence we can integrate over x to give 



\ j d* k (uE^imi^^^-^ ~ ^ ( ^vff |a ) » -<^<- (155) 

Therefore, when — 7r < < 7r the two integrals have leading terms that are equal but 
opposite, and c cont)n (t<f,) decreases faster than n~ 2 as n increases. However, we now let 
approach n. This forces B to approach vr + /4, and the integration contour moves across the 
line of poles, to its lower edge, resulting in a rapid change in the integral by an amount equal 
to the contribution of the poles. To determine the value of the new integral we increase B 
further, holding fixed at 7r + /4, so that B moves into the range ir/2 < 2B < 37r/2 where 
Re 7 is large and negative and e -271 " 7 is exponentially large; it follows that the new integral 
is negligible. Hence the integral on the left side of Eq. ( |155|) is negligible when lies on the 
second sheet. It follows that while c con t, n (^) vanishes through order 1/n 2 when t$ lies on 
the first sheet, 

, . i /3vr 2 Z/i|^(0)| 2 \ , , . 

C C ont,n(^) ~ I 25/2q I , 7T < < 37T, (156) 

when tcf, lies on the second sheet. We have now established Eqs. ( [76|) and Eqs. (|77|) . 
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